Abstract. We present the analysis of a deep colour- 
magnitude diagram (CMD) of NGC 1831, a rich star 
cluster in the LMC. The data were obtained with 
HST/WFPC2 in the F555W (~ V) and F814W (~ I) fil- 
ters, reaching 771555 ~ 25. We discuss and apply a method 
of correcting the CMD for sampling incompleteness and 
field star contamination. Efficient use of the CMD data 
was made by means of direct comparisons of the ob- 

■ served to model CMDs. The model CMDs are built by 
an algorithm that generates artificial stars from a single 

■ stellar population, characterized by an age, a metallic- 
I ity, a distance, a reddening value, a present day mass 

function and a fraction of unresolved binaries. Photo- 
metric uncertainties are empirically determined from the 
data and incorporated into the models as well. Statisti- 
cal techniques are presented and applied as an objective 
7— I ■ method to assess the compatibility between the model 
k>" I and data CMDs. By modelling the CMD of the central 
^ ' region in NGC 1831 we infer a metallicity Z = 0.012, 
^ ; 8.75 < log(r) < 8.80, 18.54 < (m - M) < 18.68 and 
'■ 0.00 < E(B - V) < 0.03. For the position dependent 
O ' PDMF slope (a = -dlog$(M)/dlogM), we clearly ob- 
serve the effect of mass segregation in the system: for 
projected distances R < 30 arcsec, a ~ 1.7, whereas 
^ '■ 2.2 < a < 2.5 in the outer regions of NGC 1831. 
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1. Introduction 

The Large Magellanic Cloud (LMC) presents three essen- 
tial characteristics that make it an excellent complemen- 
tary laboratory for studying the formation and evolution 
of galaxies and stellar systems in general: a) it is close 
to the Galaxy; b) it has markedly different morphologi- 
cal, chemical and kinematical properties from our Milky 
Way; c) it presents a large variety of stellar clusters, dis- 
playing distinct physical characteristics among themselves 
and when compared to those in the Galaxy (Wcsterlund 
1990). Due to the diversity in ages and metallicities, LMC 
clusters are found at distinct evolutionary stages (Wester- 
lund 1990, Olszewski et al. 1991). The determination of a 
cluster's present physical properties, such as density pro- 
file, shape, internal velocity distribution and its position 
dependent Present Day Mass Function (PDMF), provide 
us with essential links needed to assess the role of grav- 
itational dynamics, including effects of mass segregation 
and stellar evaporation (Hcggie & Aarseth 1992, Spurzem 
& Aarseth 1996, de Oliveira et al. 2000). Thus, once these 
present properties are known, modelling techniques like N- 
body simulations allow us to recover the initial conditions 
under which clusters formed (Goodwin 1997, Vesperini & 
Heggie 1997, Kroupa et al. 2001). In this context, the ini- 
tial mass function (IMF) and its possible universality, are 
key pieces in the study of stellar contents of distant galax- 
ies (Kroupa 2001). 

Through the analysis of colour-magnitude diagrams 
(CMDs) one has a great variety of physical information 
about a cluster. Isochrone fits help constraining the sys- 
tem's age, metallicity, distance and reddening. Further- 
more, observational luminosity functions (LFs) have al- 
lowed derivation of the stellar mass function (Elson et al. 
1995, De Marchi & Paresce 1995, Santiago et al. 1996, 
Piotto et al. 1997, de Grijs et al. 2002a). However, the 
transformation of luminosity into mass depends on age 
and metallicity, the uncertainties in these parameters be- 
ing therefore incorporated into the inferred mass func- 
tion. Besides, the theoretical uncertainties in the mass- 
luminosity relation itself, specially for the low-mass stars 
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(Piotto et al. 1997, Baraffe et al. 1998), added to the ef- 
fect of unresolved binarism, further hampers real mass 
function determination through observational luminosity 
functions. 

From both observational and theoretical points of 
view, the advent of the Hubble Space Telescope (HST), 
associated with the constant improvement in theories of 
stellar interiors, atmospheres and evolution, require ever 
more sophisticated methods of CMD analysis. In this con- 
text, computational modelling techniques have allowed a 
wider use of CMDs as tools to constrain physical prop- 
erties of stellar populations and systems. Model CMDs, 
along with statistical methods of comparing them to ob- 
served ones, have been useful means to investigate the 
star formation history within a galaxy (Gallart et al. 1996, 
Gallart et al. 1999, Hernandez et al. 1999, Holtzman et al. 
1999, Hernandez et al. 2000) or to constrain structural pa- 
rameters and the stellar luminosity function in the Milky 
Way (Kerber et al. 2001). 

With these issues in mind, our work aims at extract- 
ing as much physical information as possible about rich 
LMC clusters, by means of objective comparison of their 
observed CMDs with artificial ones. The idea is to simul- 
taneously infer PDMF slope, age, metallicity, distance, 
reddening and unresolved binary fraction for each system 
studied. The present work introduces the techniques devel- 
oped for that purpose and shows the results for NGC 1831, 
one of the richest LMC clusters for which we have deep 
HST data. 

Previous works are evidence of the large difficul- 
ties in extracting physical parameters for NGC 1831. 
Techniques based on CMDs from CCD photometry 
(Mateo 1987,1988; Chiosi 1989; Vallenari et al. 1992; 
Corsi et al. 1994), integrated spectroscopy or colours 
(Bica et al. 1986; Meurer et al. 1990; Cowley & 
Hartwick 1992; Girardi et al. 1995) and spectroscopy 
of individual giant stars (Olszewski et al. 1991) were 
employed with this aim, constraining the values of the 
main parameters: 8.50 < log(r) < 8.80 (300 < r < 650 
Myr); 0.002 < Z < 0.020 (-1.00 < [Fe/H] < 0.01); 
0.00 < E(B - V) < 0.07. For the distance to NGC 1831, 
there are not reliable determinations, the standard 
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procedure being the adoption of typical values of the 
intrinsic distance modulus, (m — M)o, for the LMC 
centre. In this aspect, the most reliable estimate seems 
to be that of Panagia et al. 1991, (m - M) = 18.51, 
since it is based on purely geometrical arguments applied 
to high quality imaging and spectral data on supernova 
SN1987a. In terms of dynamical evolution for this system, 
Elson et al. (1987) estimated 6.5 ^ log(£ cross ) ^ 7.0 and 
9.6 log(t r h) ~ 10.0 for the crossing time and two-body 
relaxation time, respectively, the range quoted being 
due to different mass-luminosity relations. Comparing 
with its estimated age values, these results suggest that 
NGC 1831 is a system dynamically well mixed, but 
not totally relaxed through stellar encounters. Hence, 
NGC 1831 is sufficiently old to have suffered mass segre- 
gation, affecting the PDMF slope at different distances 
to its centre, but perhaps still young enough that the 
initial conditions could be preserved in its outer regions. 
Similarly, external effects may not have had enough time 
to affect the cluster dynamics either. 

One of the main objectives of this paper is to verify 
the effect of mass segregation in NGC 1831, quantifying 
the variation in PDMF slope with projected distance to 
the cluster's centre. This determination may yield strong 
links to IMF reconstruction efforts based on N-body mod- 
els. The paper is outlined as follows: in Sect. 2 we present 
the data and the methods of accounting for sample in- 
completeness and field star contamination in the observed 
CMD; in Sect. 3 we present the algorithm used for CMD 
modelling and the statistical tools used for model vs. data 
comparisons; in Sect. 4 we discuss control experiments 
used for to verify the validity of the method; finally, in 
Sect. 5 we present the results for the NGC 1831 data, 
which are discussed in Sect. 6. 

2. The observed CMD 

We have data taken with the Wide Field and Planetary 
Camera 2 (WFPC2) on board HST for 8 rich LMC clus- 
ters and nearby fields. These data are part of the GO7307 
project, entitled "Formation and Evolution of Rich Star 
Clusters in the LMC" (Beaulieu et al. 1999, Beaulieu et 
al. 2001, Johnson et al. 2001). For each cluster, images 
were obtained using the F555W (V) and F814W (I) fil- 
ters. Most of the photometry had been previously carried 
out: cluster stellar LFs were built and analyzed by Santi- 
ago et al. (2001), de Grijs et al. (2002b); field stellar pop- 
ulations were studied by Castro et al. (2001). Exposure 
times, field coordinates, image reduction and photometry 
processes are described in detail by those authors. 

Fig. |l| shows the observed WFPC2 CMD for stars in 
the direction of NGC 1831 in panel (a) (hereafter the on- 
cluster field) and for a nearby field (hereafter the off- 
cluster field) studied by Castro et al. (2001) in panel (b). 
The on-cluster sample shown here is the final composition 
of the CEN and HALF images described by Santiago et 




Fig. 1. The on-cluster (a) and off-cluster (b) WFPC2 
CMDs for NGC 1831. The former contains 7801 stars ob- 
served in the cluster's direction while the latter has 2030 
stars located in a field 7.3 arcmin away from the cluster's 
centre. 



al. (2001). These have the Planetary Camera (PC) centred 
on the cluster's centre and half-light radius, respectively. 
The off-cluster field is located at about 7.3 arcmin away 
from the cluster's centre. A clear cluster main sequence 
(MS) is visible in the figure, stretching from TO555 ~ 18.5 
down to TO555 ~ 25. The cluster MS turn-off is also clearly 
visible at the upper MS end. Notice, however, that satura- 
tion becomes a problem in the HALF field for m 555 < 19 
(TO555 ~ 17.8 for the CEN field). Hence, all our subsequent 
analysis will be based on the CMD fainter than this limit. 
A branch of evolved stars is seen as well, especially in the 
range 18 ^ TO555 ~ 19, where the cluster red clump is lo- 
cated. The subgiant branch at fainter magnitudes is due to 
field contamination and is largely made up of older (r > 3 
Gyr) stars. 

The on-cluster data suffer from two important effects: 
sample incompleteness and contamination by field stars. 
Our CMD modelling algorithm does not incorporate such 
effects. Therefore it is crucial to adequately correct the ob- 
served CMD for them in order to place models and data 
on equal footing. Quantifying systematic and random pho- 
tometric uncertainties and either correcting for them or 
applying them to model CMDs is extremely important as 
well, as they are responsible for most of the observed CMD 
spread. These data corrections are the subject of the next 
subsections. 
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2.1. Random photometric uncertainties 

The random photometric uncertainties are the main 
source of spread in our HST/WFPC2 CMDs. Therefore, 
a suitable assessment of these uncertainties in both filters 
is crucial for correctly incorporating this effect into the 
artificial CMDs. 

For the on-cluster field, two independent photometric 
measurements were available for a fraction of the stars 
due to the overlap region imaged by both the HALF and 
CEN fields (Santiago et al. 2001). Thus, we used the stars 
belonging to this overlap region to estimate the typical 
photometric uncertainties in the data. For each filter and 
at each magnitude bin, we estimated the dispersion, a', 
in the distribution of differences between the independent 
magnitude measurements. For simplicity we assume that 
cr' is the composition of two equally-sized realizations of 
photometric error, a. Thus, we get a' 2 — 2a 2 , and there- 
fore 
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We emphasize that the two filters were treated as ab- 
solutely independent. As a consequence, the uncertainty 
in colour 771555 — msn, c co i our , will be the quadratic sum 
of the individual uncertainties: 



Fig. 2. Estimated photometric uncertainties. Panel (a) 
shows the 0555 (solid circles) and <7si4 (crosses) values as 
a function of apparent magnitude. Panel (b) shows the 
^colour calculated for the MS stars as a function of m.555. 



"colour — "555 °814- 

Fig. H shows the derived uncertainties for MS stars in 
the two filters, 0555 and <78i4j an d for the colour, a co \ om , 
using the expressions above. 

For the off-clusters stars, we did not have two inde- 
pendent and overlapping WFPC2 images and, as a re- 
sult, we could not apply the same method to quantify 
their photometric uncertainties. The solution found was 
to employ the same estimate as for the on-cluster stars. 
As the off-cluster image is deeper and sparser than the 
on-cluster ones, we can expect that this approximation 
yields an overestimate of the photometric uncertainties in 
the off-cluster data. However, the off-cluster CMD is used 
only for statistical subtraction of field contamination from 
the on-cluster CMD. We will see later that this particular 
correction technique is fairly independent of the estimated 
photometric uncertainties. 

2.2. Systematic photometric uncertainties 

Systematic effects in WFPC2 data have been found by 
several authors. Johnson et al. (2001) measured an ex- 
posure time effect varying from 0.01 to 0.06 for differ- 
ent WFPC2 chips and filters in images of NGC 1805 and 
NGC 1818, as part of this project, de Grijs et al. (2002b) 
finds similar trends, but of slightly larger amplitude for 
the same data. Previously, Casertano & Muchtler (1998) 
found an exposure time effect but in the opposite sense 
as Johnson et al. (2001). Colour shifts of ~ 0.04 have also 



been measured between different WFPC2 chips, possibly 
due to errors in CTE and aperture corrections or in zero- 
points (see also Johnson et al. 1999). 

Any photometric biases, as a function of either expo- 
sure time or chip, must be eliminated from our observed 
CMD, since the model CMDs which will be compared to 
it do not incorporate them. 

As there does not seem to exist a consensus on the 
corrections to be applied, our approach was to empirically 
measure such biases and to apply appropriate shifts to 
the data when necessary. We searched for both exposure 
time and chip vs. chip effects. No systematic effects were 
found in the off cluster data. The main source of bias 
in the on-cluster data was found to be an offset between 
the PC and the Wide Field Camera (WFC) chips in the 
sense that MS stars with the same 771555 magnitude tend 
to be bluer by 0.05-0.10 mag when imaged with the PC 
than with the WFC. This applies to both HALF and CEN 
images. As the PC in the HALF image is centered on the 
cluster half-light radius, this effect is unlikely to be due to 
differences in crowding. In order to correct the data for this 
effect we first defined MS fiducial lines, taking the median 
777555 — 777814 colour at different magnitude bins. This was 
done separately for each chip and each image (CEN or 
HALF). We noticed that the WFC MS lines were more 
stable, always occupying loci in the CMD very close to 
each other. Thus, we transposed the PC fiducial lines to 
the corresponding WFC locus. The corrections are shown 
in Fig. pi Panel (a) shows the uncorrected PC and WFC 
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according to the position of each observed star and its 
previously computed completeness weight. In order to fill 
the observed CMD with the missing stars, we first de- 
fined a fiducial line for the MS. As described in Sect. 2.2, 
this line was defined by taking the median 771555 — m-su 
colour at different magnitude bins. Given the i th MS star, 
its position along the fiducial line is provided by its 777.555 
magnitude and the corresponding 777555 ~~ "7,814 colour. If 
its completeness weight is Wi, then Wi — 1 extra stars 
were spread out from its position along the MS line tak- 
ing the measured photometric uncertainties (as described 
in Sect. 2.1) into account. We assumed a Gaussian dis- 
tribution of uncertainties both in 777555 magnitude and 
W555 — 777814 colour. Sample completeness falls to less than 
50% for 777555 ~ 23.8. We completed the CMD down to 
W555 = 24.5 and then cut it at 777555 = 23.5, therefore 
avoiding boundary effects. A complete sample of stars 
both in number and in position resulted from this method. 

As for the off-cluster CMD, it is complete at least down 
to 777555 = 24.5 (Castro et al. 2001). 



Fig. 3. MS fiducial lines for each chip (PC or WFC) and 
image (CEN or HALF) as indicated: (a) original MS lines; 
(b) MS lines after applying the correction discussed in the 
text. 



(a mean locus of the 3 chips) lines for both the CEN and 
HALF images. Panel (b) shows the corrected fiducial lines. 

2.3. Selection effects 

2.3.1. Sample incompleteness 

Sample incompleteness occurs essentially due to two fac- 
tors: overlapping of stellar profiles due to crowding and 
background noise. Therefore, in a cluster incompleteness 
will depend not only on apparent magnitude, but also on 
spatial position. Faint stars in dense stellar regions are 
the ones that suffer most from this effect. Completeness 
effects in the on-cluster data were previously discussed 
and measured by Santiago et al. (2001). These authors 
carried out experiments where artificial stars were added 
to the cluster image and subjected to the same sample 
selection as the real stars. The completeness c of each 
real star was estimated as the fraction of artificial stars 
of similar magnitude and location which were successfully 
recovered in the experiments. The estimated weight w for 
each star is then simply given by the inverse of its com- 
pleteness (w = 1/c). This weight corresponds to the total 
number of stars similar to the observed one which should 
be detected and measured in an ideal image. 

Assigning a weight to each observed star is enough for 
the sake of luminosity and mass functions. However, cor- 
recting a CMD for incompleteness requires an extra step, 
namely to place an extra number of stars on the CMD 



2.3.2. Field star subtraction 

Field star subtraction from the cluster sample is carried 
out in two separate steps. We first cut-off all stars in the 
on-cluster CMD which are much farther from the MS than 
expected by photometric errors. We therefore eliminate 
all evolved stars, as well as objects which are likely to be 
foreground stars or remaining non-stellar sources in the 
sample (distant galaxies, spurious image features, etc). 

The second step involves the statistical removal of field 
stars located along the MS. We thus compare the dis- 
tribution of stars in the on-cluster CMD to that of the 
off-cluster CMD. The comparison method is based on the 
hypothesis that the positions of the off-cluster stars repre- 
sent the most likely positions for field stars on any similar 
CMD. We thus try to estimate the probability of each 
on-cluster star to be a field star and according to this 
probability we randomly remove stars from the on-cluster 
CMD. 

We consider pairs of on-cluster/off-cluster stars. For 
each pair we compute the expected scatter in CMD po- 
sition between the pair members under the assumption 
that they are independent photometric realizations of the 
off-cluster star. We will then have 

°c,555 = \Z2cToff_555 

and 

0~c, colour V^2<Toff, colour 

respectively for the expected scatter in 777555 magni- 
tude and 777555 — »^8i4 colour, where the "off" subscript in 
the expressions above refers to the off-cluster pair mem- 
ber. 
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For the i th off-cluster CMD star we then consider the 
Ni on-cluster stars inside a 3<t Cj 555 x 3cr c , co i our box cen- 
tered on it. Using a Gaussian error distribution in mag- 
nitude and colour, we estimate the probability pij that 
the j th on-cluster CMD star, inside this box, is a second 
photometric measurement of the i th field star. Therefore, 



-(magi - magj) 2 
Ptj oc exp[ — — \exp 



-{colour i — colour jY 



2(c C! 555) 2 * 2((T C!CO i our ) 2 

where the normalization of pij is such that 

J2pi,j = 1 

3=1 

Doing the same for all N a s off-cluster CMD stars, we 
estimate the probability Pj that the j th on-cluster CMD 
star is any one of the N Q g field stars. Hence, 



P, = Vftj, 



4=1 



where 



3=1 

and iV on is the total number of stars in the on-cluster 
CMD. In practice, not all off-cluster stars will have on- 
cluster stars within their 3cr Cj 555 x 3cr c co i our boxes. The 
actual off-cluster stars taken into account will therefore 
be N^ B < N oS . 

Based on the Pj probabilities, and scaling the number 
of field stars to the solid angle of the on-cluster field, we 
randomly remove 



N 6t 
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N' ° a 

^off 
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stars from the on-cluster CMD, where Q on and Q s are 
the solid angles covered by the on-cluster and off-cluster 
fields, respectively. 

Fig. [| shows the results of correcting an observed 
NGC 1831 CMD for incompleteness and field contamina- 
tion. Panel (a) shows the original CMD obtained from 
photometry, excluded of non-MS stars and corrected only 
for the systematic photometric effects described in Sect. 
2.2. Panel (b) presents the complete CMD, i.e., corrected 
for sample incompleteness and cut at 771555 = 23.5; the 
extra stars added by completeness correction represent 
about 23% of the total within this magnitude range. Panel 
(c) shows the 89 stars (~ 3%) in the on-cluster CMD that 
were considered as LMC field stars in the field star sub- 
traction process; finally the clean and final cluster CMD 
is shown in panel (d). 

We tested the field star subtraction algorithm for dif- 
ferent box sizes and assumptions regarding the photomet- 
ric scatter. The results are insensitive to the details in the 
algorithm. 



rn — 1 — 1 — 1 — 1 — r 

Complete CMD 
8737 stars 




Fig. 4. Correcting the on-cluster CMD for incompleteness 
and field stars contamination: (a) the observed CMD along 
the MS; (b) the CMD corrected for sample incompleteness 
and cut at 771555=23.5. The MS fiducial line for the cluster 
is shown in white ; (c) on-cluster stars considered as LMC 
field stars according to the algorithm described in the text; 
(d) final CMD with only stars presumed to belong to the 
cluster. 



3. CMD modelling and statistical tools 

3.1. CMD Modelling 

We model the MS of a cluster as a single stellar popula- 
tion, characterized by stars of the same age and chemical 
composition. The first step is the choice of an isochrone, 
which defines a sequence of magnitude and colour as 
a function of mass for this population. For the present 
work we used Padova isochrones (Girardi et al. 2000) be- 
cause they present masses inside the observed MS mass 
range [M ;$ 2.5 M@) and are expressed in the vegamag 
WFPC2 photometric system. Notice that Padova mod- 
els assume convective overshooting in the stellar interi- 
ors. This assumption may influence age determinations, 
specially when based in the position of turn-off and He- 
burning stars of younger populations (Testa et al. 1999, 
Barmina et al. 2002). However, as NGC 1831 is at least 
several Myr old and our modelling makes use of the entire 
MS, we believe that these model uncertainties tend to be 
of smaller importance in our analysis. 

The free model parameters are: metallicity (Z), age 
(log(r)), intrinsic distance modulus (777 — M)q, reddening 
(E(B-V)), PDMF slope (a = — d log $(M) /d log M) and 
unresolved binary fraction (/bin)- The PDMF was consid- 
ered as a power-law ($(M) = $ M _Q! = diV/dM), where 
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the only free parameter is the slope a. /bin is defined as 
the systemic binary fraction, / bin = N hin /(N hia + N sing ), 
where _/Vbi n and N S i ng are, respectively, the number of un- 
resolved pairs and single stars. 

The process of generating artificial stars works as fol- 
lows: 

(1) we fix Z and log(r) for the stellar population by 
means of a chosen isochrone; 

(2) we randomly draw a stellar mass according to the 
PDMF and get the absolute magnitudes in the two desired 
filters through the mass-luminosity relation given by the 
isochrone; 

(3) for /bin randomly chosen cases, we repeat step (2), 
representing a companion star in a binary system, and 
combine the two luminosities in both filters; 

(4) we apply the intrinsic distance modulus (m — M)q 
and reddening vector( J 4y J E(B — V)) to the system, defin- 
ing its theoretical CMD position. For this purpose, we use 
Ry = Ay/E(B — V) = 3.1 and the photometric trans- 
formation to the vegamag WFPC2 system according to 
Holtzman et al. 1995a; 

(5) we introduce the photometric uncertainties by 
spreading the star with a Gaussian distribution of errors 
with 0555 and asu as empirically determined (see Sect. 
2.1). This yields observational versions of the magnitude 
and colour; 

(6) finally, we verify if the star is inside the MS obser- 
vational ranges in magnitude and colour defined for the 
data and throw it away if it is not. 

For each model realization we generate the same num- 
ber of MS stars as observed in the real CMD, corrected 
for sample incompleteness and field contamination, and 
located inside the 19.0 < 771555 < 23.5 range. This range 
in apparent magnitude corresponds to 0.5 M555 5.0 
and 0.9 < M < 2.3 M . 

The best models are chosen by a direct comparison of 
the observed CMD with the artificial ones. The statistical 
tools used in this comparison are presented in the next 
section. The model vs. observation comparison strategy is 
as follows: 

(1) the global parameters for the cluster, log(r), Z, 
E(B — V) and (m — M)o, are determined using the CMD 
of the central cluster region (R < 30 arcsec, where R is the 
projected distance from the cluster's centre), where field 
contamination and statistical noise are minimized (see Ta- 
ble0); 

(2) for the best combinations of the global parameters, 
the position dependent parameters a and /bin are then 
derived using the CMDs in concentric rings of variable 
radii. 

Table ^ shows important information about the cluster 
regions used in this modelling process. The first column 
gives the range in R, Col. (2) the original number of CMD 



Table 1. Number of points in the different cluster regions 
whose CMDs are used in the CMD modelling as described 
in Sect. 3.1. 



Region (arcsec) 


observed 


complete 


field stars 


cluster 


< R < 30 


2221 


2737 


89 


2648 


< R < 15 


1220 


1506 


27 


1479 


15 < R < 30 


1001 


1216 


62 


1154 


30 < R < 60 


1692 


1972 


240 


1732 


R > 60 


1673 


1780 


765 


1015 



stars in the 19.0 < 771555 < 23.5 range, Col. (3) the num- 
ber of stars in the completeness corrected CMD (see Sect. 

2.3.1) , Col. (4) the number of assigned field stars (Sect. 

2.3.2) and, finally, Col. (5) the number of stars in the CMD 
used in the modelling process. Notice that, not unexpect- 
edly, field contamination becomes a serious issue for the 
outermost region, since statistical removal of field stars 
reduces the CMD numbers by about 40%. On the other 
hand, crowding in the central regions yields larger photo- 
metric incompleteness, as reflected by the clear increase 
in the number of stars between Cols. (2) and (3). 

This modelling strategy allows an efficient and system- 
atic search of best fit models in a 6-dimensional parame- 
ter space and makes use of the 2-dimensional information 
provided by the CMD data. Furthermore, this strategy 
naturally splits the parameters into those that define the 
position of the MS in the CMD plane (the global ones) and 
those that influence the way stars are distributed within 
the MS locus (the position dependent ones). 

Fig. U shows four CMDs for the central cluster re- 
gion. The one in the bottom right (panel d) is the real 
data, whereas the other three are realizations from differ- 
ent models generated by the modelling process described 
above. The input model parameters are shown in each 
panel. The three model CMDs do in fact look different, 
their MS ridge lines having different shapes and occupy- 
ing different positions along the CMD plane. 

3.2. Statistical tools 

One of the main goals of this paper is to establish an 
objective comparison method between models and data. 
This required developing and applying statistical criteria 
that discriminate the model CMDs that most adequately 
reproduce the observed one. Ideally these comparison cri- 
teria should be both simple and easy to implement but 
yet make use of as much information provided by the bidi- 
mcnsional colour-magnitude plane as possible. We stress 
that these methods, in principle, arc not restricted only 
to CMD analyses, but may be applied to the comparison 
of any two bidimensional distributions of points. In simi- 
lar context as in this work, statistical techniques of CMDs 
comparison have been successfully applied by Kerber et 
al. (2001), Lastennet & Valls-Gabaud (1999), Saha (1998), 
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Fig. 5. CMDs for the cluster central region: panels a, b 
and c show realizations from the modelling process while 
panel d shows the result of the real data treatment (Sect. 
2). The input model parameters are indicated in each 
panel, where a = 2.30 and /b; n = 0.30 for all models. 



Hernandez et al. (1999), allowing a model reconstruction 
of the main CMD features displayed by the component 
stars within a galaxy or a cluster. 

Three statistics of CMD comparison were considered 
in this paper: S 2 , PSS and L. The first two are explained 
in more detail by Kerber et al. (2001). S 2 is essentially 
a dispersion between model and data points in the CMD 
plane; PSS is proportional to the joint probability that 
the two CMDs being compared are drawn from the same 
population. 

As for the L statistics, it is an empirical version of the 
likelihood statistics described and used by Hernandez et 
al. (1999). For each model, 300 realizations with the same 
number of artificial stars as the real data (hereafter N Q b s ) 
were generated. Dividing the CMD plane into N\> boxes, 
the model probability of one star, randomly chosen from 
any of these 300 realizations, to be inside the k th box is 
given by pk — A^/(300 iV bs)- Afc is the sum of stars from 
all realizations which fall in the k th box. Thus, pk may 
be interpreted as a probability function along the CMD 
plane. 

The likelihood L of a given model is then defined as 

N obB 

L = JJ Pobs,i 
i=l 

or 



N obs 

l0 § L = l0 g(Pobs,i), 
i=l 

where the product and sum are over the -/V bs observed 
stars, and p bs,i is the model probability function evalu- 
ated at the CMD position of the i th observed star. 

So, for each model we have three distinct statistical 
values. In order to establish the best models we build diag- 
nostic diagrams (hereafter DDs), which are planes where 
we confront any two of these statistics. The best models 
will naturally have large L and PSS, and small S 2 values. 
The method was tested by means of control experiments, 
that are shown in Sect. 4. 

3.3. The model grids 

The model input values for log(r), Z, (m — M) , 
E(B — V) were chosen in order to bracket those 
found in the literature. In this respect, the web page 
www.ast.cam.ac.uk/STELLARPOPS/LMCdatabase by 
Richard de Grijs was very useful as it includes a very 
large compilation of parameter values and references on 
the 8 LMC rich clusters imaged by the GO7307 project. 

In accordance with this database, we set the range of 
possible values for each physical parameter and created 
a regular model grid within this range. We expect this 
systematic grid to prevent biases in the parameter values 
determination. 

Using the cluster central region (0 < R < 30 arcsec), 
we explored the following space defined by the global pa- 
rameters: 

Z = 0.004,0.008,0.012 

log(r/yr) = 8.70,8.75,8.80 

(m - M) = 18.30, 18.40, 18.50, 18.60, 18.70 

E(B - V) = 0.00,0.02,0.04,0.06,0.08,0.10 

The position dependent parameters were kept fixed as 
a = 2.30 and /b; n = 0.30. Therefore this initial grid has 
270 models. 

Once the CMD comparison statistics defined in Sect. 
3.2 are computed for each model, the DDs are built and 
the best models are identified, we investigate the dynami- 
cally affected, position dependent parameters by generat- 
ing artificial CMDs to be compared to the observed CMDs 
within the concentric rings. In this second step we explored 
the following bidimensional parameter space: 

a = 1.40,1.70,2.00,2.30,2.60 

/bin = 0.25,0.50,0.75 

A total of 180 models were built for each ring, the only 
difference between one ring and another being the number 
of artificial stars. 
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4. Control Experiments 

As mentioned in Sect. 3.2, we tested the validity of our sta- 
tistical methods using control experiments. These experi- 
ments consist of drawing one realization of some specific 
model and calling it the "observed CMD" . We then verify 
if the DDs recover the generating model (hereafter input 
model) as the best one describing the "observed CMD" . 

Figs. ^ and show the results of a control experi- 
ment involving the 270 models to be latter applied to the 
NGC 1831 central region. All panels show DDs of logS 2 
vs. log L, each point in the DD representing a particu- 
lar model. The panels on the right are blow-ups of those 
on the left, showing in detail the region where the best 
models are located; this region corresponds to larger log L 
and smaller log S 2 values. The different symbols in a sin- 
gle panel are coded according to the values of one of the 
four global model parameters (Z, log(r), (to — M)q and 
E(B — V)), therefore allowing the effect of varying each 
parameter to be separately assessed. Notice that the fig- 
ures do not show the entire model grid in order to avoid 
cluttering. The grid regions discarded from the DDs are 
those of systematically high logS* 2 and low logi values. 
The "observed CMD" was taken to be a realization of the 
model with Z = 0.008, log(r) = 8.75, (m - M) = 18.50 
and E(B — V) = 0.06. This input model is shown as the 
large star in the blowup panels. 

A tight correlation between the two statistics is clearly 
seen in all panels, adding reliability and stability to the 
choice of the best fitting models. The control experiments 
also show that one is in fact capable of recovering the input 
model based on the values of the statistics used: it is the 
model with largest logL and smallest logS" 2 , as ideally we 
would expect. 

Another important result from Figs. ^ and [7] is the 
combining and/or canceling effect of some global parame- 
ters, yielding models of comparable quality. As an exam- 
ple, the effects of metallicity Z and reddening E(B — V) 
tend to cancel each other. Some models with Z (E(B~V)) 
lower (higher) than the input value, along with some high 
Z and low E(B — V) ones, rank among the best models 
in the DDs. This degeneracy in the DDs is not surprising 
since the effect of increasing Z is to make stars redder and 
fainter at a given mass, roughly opposite to the effect of 
decreasing E(B — V). 

Fig. H presents similar DDs as in Figs. | and 0, but 
shows the model grid to be applied to the concentric re- 
gions of NGC 1831 (in a total of 180 models). The sym- 
bols now indicate different values of the PDMF slope a 
(panels (a) and (b)) and binary fraction /bin (panels (c) 
and (d)). As before, the panels on the right are blow- 
ups of the ones on the left, showing the best models only. 
The input model (large star), in this case, is the one with 
a = 2.00, /bin = 0.50 (and Z = 0.012, log(r) = 8.75, 
(m - M) = 18.60, E(B — V) = 0.02). It is again located 
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Fig. 6. DDs resulting from the control experiment for 
the cluster's central region, showing the effects of vary- 
ing metallicity (panels (a) and (b)) and age (panels (c) 
and (d)). The symbols are as indicated in the panels on 
the left. The panels on the right show the best models in 
detail and use the same symbol notation. 
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Fig. 7. DDs resulting from the control experiment for the 
cluster's central region, showing the effects of varying dis- 
tance (panels (a) and (b)) and reddening (panels (c) and 
(d)). The symbols are as indicated in the panels on the 
left. The panels on the right show the best models in de- 
tail. 
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Fig. 8. DDs resulting from the control experiment for the 
15 < R < 30 arcsec region. Panels (a) and (b) (the latter is 
a blowup of the former) show the effect of varying PDMF 
slope a, whereas panels (c) and (d) (the latter is a blowup 
of the former) show the effect of varying /bin- The input 
model is represented by the large star. 

at the extreme upper left corner of the DD, confirming the 
applicability of our statistical approach. 

However, the panels on the right show that the three 
best models have the same a (= 2.00) but different /bin 
values, revealing a larger difficulty in determining the lat- 
ter than the former. This occurs because the effect caused 
by binaries is of the same order as or smaller than the 
photometric uncertainties in our WFPC2 CMDs. Conse- 
quently, as will be discussed in Sect. 5.2, the /bin determi- 
nation by means of our CMDs becomes a difficult task. 

Notice that, in comparison with the DDs for the cen- 
tral region, the DDs in Fig. ||| present a larger dispersion. 
This is caused by the much smaller number of stars used 
in the set of models with varying a and /bi n ; in this case, 
all model CMDs (including the "observed" one) have 1154 
stars, therefore mimicking the situation of the second con- 
centric region to be studied in NGC 1831 (see Table [l]). 

Fig. [9] shows DDs involving the PSS statistics. The 
upper panels show PSS vs. log L and PSS vs. log S 2 plots 
for the central region. The lower panels show the same 
plots for one of the concentric rings. The same models 
as in Figs. [| 0, and || are depicted but without symbol 
coding as a function of parameter values. The correlation 
between PSS and the other statistics is again quite tight. 
In fact, the results based on the DDs are insensitive to the 
particular choice of statistics to be plotted. This is a very 
important result, since it further enhances the reliability 



Fig. 9. DDs involving the PSS statistics in control exper- 
iments. The upper (lower) panels show PSS vs. log L and 
PSS vs. logS* 2 plots for the region inside < R < 30 
(15 < R < 30) arcsec. 

of our statistical CMD modelling techniques. For the sake 
of simplicity, we hereafter restrict ourselves to log L vs. 
log S 2 DDs only. 

5. Results 

5.1. Central region 

In this section we model the CMD of stars belonging to the 
region within R < 30 arcsec from the centre of NGC 1831. 
In accordance with Santiago et al. (2001), this corresponds 
roughly to R ,$ 7 pc or within 2 half-light radii. As indi- 
cated in Table [l], there are 2221 stars in this region in 
the magnitude range 19.0 < TO555 < 23.5. As also men- 
tioned earlier, this central cluster region suffers from in- 
significant contamination by LMC field stars. On the other 
hand, due the high stellar density, incompleteness effects 
become more important. The chosen faint magnitude cut- 
off represents the magnitude at which completeness falls 
at 50%, in an attempt to reduce the relevance of this effect 
on our results. 

The DDs for the model vs. real CMD comparison are 
presented in Figs. [l(] and [ll]. As in the control experi- 
ments, the correlation between log L and log S 2 is quite 
noticeable, the best models being again at the upper left 
region in the panels. These figures follow the same con- 
ventions and notations as Figs. [6] and [?], therefore allowing 
us to assess the effect of varying each global cluster pa- 
rameter separately. For example, panels |l0](a) and |To|(b) 
clearly show that the best models have Z = 0.012. The 
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Fig. 10. DDs resulting from the CMD of the cluster's cen- Fig. 11. DDs resulting from the CMD of the cluster's cen- 
tral region, showing the effects of varying metallicity (pan- tral region, showing the effects of varying distance (panels 
els (a) and (b)) and age (panels (c) and (d)). The symbols (a) and (b)) and reddening (panels (c) and (d)). The sym- 
are as indicated in the panels on the left. The panels on bols are as indicated in the panels on the left. The panels 



the right show the best models in detail. 



on the right show the best models in detail. 



effects of the other parameters are not as striking as in 
the case of metallicity. Yet, panels |lO|(c) and 0(d) favour 
an age in the range 8.75 < log(r) < 8.80. Likewise, the 
three lowest values of intrinsic distance modulus (two of 
which are not even shown in the figure) can essentially 
be ruled out (panels 0(a) and 0(b)), placing NGC 1831 
near or beyond the distance to the LMC centre. Finally, 
E(B — V) < 0.02 is favoured by our modelling approach 
(panels 0(c) and 0(d)). Notice that this latter parameter 
is again anti-correlated with Z, in the sense that the mod- 
els in the blowup panels with Z — 0.008 (circles in panel 
0(b)) have larger E(B — V) values (circles or crosses in 
panel 0(d)). 

These best choices of the global parameters, as dis- 
cussed in Sect. 3.1, besides being useful constraints by 
themselves, serve as model inputs to the position depen- 
dent parameters a and /bin, whose modelling is based on 
the concentric regions. 

5.2. Concentric regions 

We now analyze the NGC 1831 CMDs in the concentric re- 
gions listed in Table |l|. As just mentioned, the free model 
parameters in this case are the PDMF slope a and the un- 
resolved binary fraction /bin- For all regions we corrected 
the observed CMD to the effects of sample incomplete- 
ness and field star contamination. Fig. O shows their final 



CMDs. As mentioned before, Table [j] shows the number of 
points in all the steps along the data correction process. 

In the modelling, we apply the same model grid for 
every region, only changing the number of artificial stars 
generated. The results from the previous section were used 
to restrict the possible values of the clusters global param- 
eters, allowing us to investigate the dynamically depen- 
dent parameters for the system. 

Figs. and present the DDs resulting from model 
vs. data CMD comparison in the four regions. The differ- 
ent symbols in each panel correspond to different values of 
a. As usual, the panels on the right show the best models 
in detail. The results for the two innermost regions are 
presented in Fig. Panels (a) and (b) indicate that the 
best models for innermost region (R < 15 arcsec) have 
1.40 < a < 1.70. In the ring 15 < R < 30 arcsec, shown 
in panels (c) and (d) of the same figure, the best models 
present similar a values: 1.40 < a 1.70. 

In Fig. we have the results for the two outermost 
rings. Panels (a) and (b) show the DDs for the stars with 
30 < R < 60 arcsec. Now, there is evidence for a sharp 
change in the PDMFs slope: the best models have 2.30 < 
a < 2.60. Finally, the DD for the most peripheric ring has 
a large dispersion (panels 0(c) and 0(d) ). This spread 
is also seen in the best PDMF slope: values in the range 
2.00 ^ a < 2.60 are seen in the upper left corner of the 
DD. This may reflect large uncertainties in the field star 
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Fig. 12. The final NGC 1831 CMDs in four concentric re- 
gions: (a) < R < 15 arcsec; (b) 15 < R < 30 arcsec;(c) 
30 < R < 60 arcsec; (d) R > 60 arcsec. 

Table 2. Best estimates of the slope a and its uncertainty 
as a function of distance to the centre. 



Region (arcsec) 


Region (pc) 


a 


Pot 


< R < 15 


< R < 3.8 


1.72 


0.15 


15 < R < 30 


3.8 < R < 7.6 


1.68 


0.19 


30 < R < 60 


7.6 < 7? < 15.2 


2.45 


0.15 


R > 60 


R > 15.2 


2.19 


0.33 



subtraction, since it leads to a large reduction of CMD 
points in this ring. 

For each ring, we assign a representative a and asso- 
ciated uncertainty using the 10 models with the largest 
values of likelihood L. The final slope is taken to be the 
weighted average value among these best models, and its 
associated uncertainty is the dispersion around the aver- 
age. The weight assigned to each model was the inverse of 
the difference in logi between the observed CMD and a 
CMD from a typical model realization. This difference is a 
measure of the discrepancy between model and observed 
CMDs. Table || lists the final a values and uncertainties 
for each ring (Cols. 3 and 4). The first 2 columns in the 
table show the regions in arcsec and parsecs (assuming 
(m - M) = 18.61 for NGC 1831 as adopted in Sect. 6), 
respectively. 

As for the fraction of unresolved binaries, /bin, the re- 
sults, as anticipated, are not conclusive in any ring, given 
the CMD spread. A more precise treatment of photometric 
errors may help constrain this parameter. 



Fig. 13. DDs resulting from the CMDs of concentric re- 
gions, showing the effect of varying a. The upper (lower) 
panels show the results for stars in the region with < 
R < 15 (15 < R < 30) arcsec. The panels on the left 
show the entire DD, while those on the right show the 
best models in detail. The symbols are as indicated in the 
panels on the right. 

6. Summary and Conclusions 

In this work we analyzed a deep CMD of NGC 1831, a 
rich LMC cluster, obtained from HST/WFPC2 images in 
the F555W and F814W filters. We inferred physical pa- 
rameters such as metallicity, age, intrinsic distance modu- 
lus, reddening and PDMF slopes by comparing the cluster 
CMD with artificial ones. We presented in detail the tech- 
niques used to build the model CMDs, which use these 
parameters as model input and take into account obser- 
vational uncertainties as in the real data. The parameter 
space explored by our regular model grids bracketed the 
values found in the literature. 

The model vs. data comparison required correcting the 
observed CMD for selection effects caused by photometric 
incompleteness and CMD contamination by stars belong- 
ing to the LMC field. We also presented the statistical 
techniques used to compare the real CMD, corrected for 
the aforementioned selection effects, to the artificial ones. 
These statistical tools allowed us to discriminate the mod- 
els that best reproduce the data. They are based on simple 
and objective statistical quantities and make use of the full 
bidimcnsional distribution of points in the CMD. 

The best parameter values inferred for NGC 1831 are 
in the ranges 8.75 < log(r) < 8.80, 18.50 < (m - M) < 
18.70, 0.00 < E(B -V) < 0.02. Using the weighted aver- 
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Fig. 14. DDs resulting from the CMDs of concentric re- 
gions, showing the effect of varying a. The upper (lower) 
panels show the results for stars in the region with 30 < 
R < 60 (R > 60) arcsec. The panels on the left show the 
entire DD, while those on the right show the best models 
in detail. The symbols are as indicated in the panels on 
the right. 

age value among the 10 best models, as described in Sect. 
5.2, we have r = 588 ± 43 Myr, (m - M) = 18.61 ± 0.07, 
E(B - V) = 0.013 ± 0.015. As for the metallicity, all 10 
best models have Z = 0.012. 

As discussed in Sect. 4, there is a coupling between 
reddening and metallicity, and therefore the determina- 
tion of the former limits the values of the latter. As a 
consequence, the derived high metallicity for NGC 1831 
brings about a low reddening determination. 

The combined values of the global parameters obtained 
in this work for NGC 1831 suggest that this cluster is: a) 
metal-richer and older than in most previous estimates; b) 
placed near or beyond the distance to the LMC centre; c) 
found in a low reddening region. 

Our estimated age and metallicity for NGC 1831 are 
in perfect accordance with the age-metallicity relation for 
LMC clusters obtained by Olszewski et al. (1991). 

The best models for the global parameters were then 
used to build artificial CMDs of regions separated accord- 
ing to distance from the cluster centre and with varying 
values of the PDMF slope and fraction of unresolved bina- 
ries. For the PDMF slope a our statistical modelling shows 
that significant mass segregation exists in NGC 1831: for 
the central regions (out to 30 arcsec ~ 7.6 pc), we derive 
a ~ 1.70, whereas a £ 2.20 for the outer regions. As for 
binarism, our results were not as conclusive. One explana- 



tion is that the spread in the CMD caused by unresolved 
binaries is of similar or smaller amplitude than the empiri- 
cally derived photometric uncertainties in the data. There- 
fore, this parameter is highly sensitive to the adopted pho- 
tometric errors, rendering its realistic estimate a task for 
yet deeper CMDs or for data for which photometric un- 
certainties may be better estimated and modelled. 

In a preliminary analysis of these data, Santiago et 
al. (2001) observed the effect of mass segregation in 
NGC 1831 with a steepening of the LF slope as a func- 
tion of distance from the centre. However, no PDMF was 
derived. A global PDMF slope for NGC 1831 is presented 
by Mateo (1988), through the conversion of the luminos- 
ity function down to V ~ 23 (M > M Q ) into a PDMF. 
The resulting slope is a <~ 4.0, therefore considerably 
steeper than our position dependent ones. Global slope 
values of other clusters from Mateo (1988) are in the range 
2.5 < a < 4.6. In contrast, Elson et al. (1989) find much 
shallower slopes 0.8 £ a £ 1.8 for another sample of rich 
LMC clusters. These earlier results are based on ground- 
based observations, for which the effects of crowding are 
much more serious than in the present work. 

Our determined position dependent PDMFs constrain 
the a values within ~ 0.2. These results can be used along 
with N-body simulations in order to recover the initial 
conditions, in particular the cluster IMF. We are applying 
the techniques shown in this paper to the others LMC 
clusters imaged with HST/WFPC2 as part of the GO7307 
project in order to infer the same physical information as 
for NGC 1831. These future results can be very useful in 
investigations on the IMF universality. 
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